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Abstract. In this review we address the uncertainties implicit in evolution- 
ary synthesis model computations. After describing the general structure of 
synthesis codes, we discuss several source of uncertainties that may affect their 
results. In particular, we discuss the uncertainties arising in the computation of 
isochrones from evolutionary tracks; those related to atmosphere models; those 
that are a consequence of the incompleteness of the input ingredients; and those 
associated with the computational aspect used in synthesis codes. We also dis- 
cuss the issue of the inclusion of distributed properties in synthesis models, an 
issue that will become relevant in the next future; as a paradigm of this case, we 
illustrate the difficulties implied by the inclusion of tracks with rotation in syn- 
thesis models. Finally, we describe several examples of the statistical approach 
to population synthesis. 

We report on the failure of the fuel consumption theorem (FCT) and the 
isochrone synthesis code to produce mutually consistent results. However, we 
argue that FCT and isochrone synthesis results are reliable for application to 
real systems in the wavelength range where they coincide. 

On the constructive side, we derive several useful survival strategies to 
bypass uncertainties. We show that single stellar populations at the turn-off 
ages of the tabulated tracks can be safely compared, as they are scarcely affected 
by the interpolation scheme used to compute isochrones. Finally, we suggest to 
use derivative quantities, such as the SN-rate, as bug detectors. 

On the recommendation side, we advocate for greater transparency and 
more documentation in synthesis modeling. We also ask stellar model makers 
to think of us and include more mass values in the tracks. 



1. Introduction 

Synthesis codes are used for a variety of scopes. Their goals range from testing 
the accuracy of evolutionary tracks to explaining the physical characteristics of 
distant galaxies, to measuring the extinction of outer galaxies, or to computing 
suitable inputs to photoionization models. They are a natural link between our 
knowledge of individual stars and the properties of stellar ensembles, and as 
such they are an invaluable tool in the analysis of stellar populations. 

Given this flexibility, it is tempting to forget their limitations; yet an aware- 
ness of the uncertainties affecting the results of synthesis codes is necessary if 
one wants to take full advantage of their power. Limitations arise from different 
issues. One is specialization: although, as a generic tool, synthesis codes can 
be used to address a variety of problems, each code is specialized in a particu- 
lar niche of the general physical problem. Stretching the result of a particular 
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code out of its natural boundaries is a frequent mistake, and statements like 
"Everybody knows that the results of synthesis codes are code-dependent" are 
unfortunate outcomes of code misuse. Other limitations are a direct consequence 
of the incompleteness of the ingredients used; still others depend on the features 
of the numerical methods employed by the code. All these uncertainties must 
be known and understood whenever synthesis models are used. 

In order to assess the uncertainties in stellar population synthesis, it is 
mandatory to understand what evolutionary synthesis codes do, which their in- 
gredients and their intrinsic assumptions are, how such assumptions are managed 
by the codes, and to which situations the results of the codes can be applied. 
Following this scheme, the structure of this review is the following. We will first 
describe the general structure of any synthesis code in section SjT]In section 
we will discuss the uncertainties and limitations of the ingredients and their im- 
pact on the interpretation of the results; we include here an extensive discussion 
of how isochrones can be obtained from evolutionary tracks, and their associated 
uncertainties. In section ^]we will discuss the problems related to numerical 
methods. Section ^T]is devoted to the challenges posed by the inclusion of stel- 
lar rotation and variability in synthesis codes. In section £jfT]we will discuss the 
interpretation of results and in section ^7]we will draw the conclusions. 

Before starting, let us apologize to those who expect to find an extensive set 
of references in this review. Most of the issues addressed here are also discussed 
by other papers reporting on the results of particular evolutionary synthesis 
codes, sometimes described at length and sometimes mentioned in just a few 
sentences. Instead of citing all these papers, we have chosen to maintain the 
number of references to a minimum preferring to discuss exhaustively the issues 
related to the unce rtainties in model s. Fo r more traditional rey i ews o n synthesis 
models we refer to iMarastonl (l2003f) and IPodov Sz Prokhorovl (20(3). We also 
refer to the papers bylCharlot et al.l (jl996t l: Ide Griis et al.l (|2003l ): lAnders et all 
(2004h: lde Griis et all (|2005T t. which discuss some issues that are not addressed 
here. 



2. Overview of the problem 

The general problem addressed by synthesis codes is the computation of the 
luminosity L to t emitted by an ensemble of iVtot stars - a stellar population. From 
a theoretical point of view, this problem can be characterized in the following 
basic ways. 

If the luminosities l\ of the individual stars are known, the total luminosity 
Ltot is obtained trivially as the sum of all the l\ values: 

iVtot 

£tot = E (!) 

This circumstance is not very frequent; its most common examples are Monte 
Carlo synthetic clusters in the theoretical domain, and resolved stellar popula- 
tions in the observational one. 

In the more usual situation in which the luminosities of individual objects 
are not known, a different approach must be adopted. In this case, it can be as- 
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sumed that the emission of the ensemble is given by the sum of the contributions 
of iVd^ different classes of stars: 



1=1 



each including a relative number of stars Wi with average luminosity li. A class 
is formed by fairly homogeneous stars, a group of contiguous classes can be 
identified with a stellar evolutionary phase 1 , and the whole of classes represents 
a particular ensemble of stars. 

As will be explained in the following, the coefficients Wi can be computed 
with the results of stellar evolution and stellar atmospheres theories weighted 
by the stellar birth-rate, a function returning the number of stars of each given 
initial mass born at a given time. 

Stellar evolution theory describes the time evolution of model stars. Its re- 
sults are quantities in the theoretical space, e.g. the bolometric luminosity Lbol, 
the effective temperature T e fj, and so on. However, the results must be com- 
pared with observable quantities, such as luminosities in a given band or spectral 
indices. Hence, a transformation from the theoretical space to the observational 
one is required: this is the domain of the theory of stellar atmospheres. Such 
transformation requires using a collection of templates of the emission of the 
considered source (i.e. an atmosphere library). In mathematical terms, this 
corresponds to defining the observable quantity li in terms of quantities in the 
theoretical space, i.e. li{L\> \,T e R, ...). 

The stellar birth rate tells us how many stars with a give n initial mass are 
born in the cluster at a given time. It is often assumed (e.g., iTinslev &: Gunr] 
1976) that the stellar birth-rate can be separated into two functions independent 
of each other: the Star Formation History (SFH), which gives the number of 
stars born in a given time, and the Initial Mass Function (IMF), which gives the 
relative number of stars born as a function of the mass. This assumption, which 
will be discusse d briefly in H2.3.1 is the strongest ad hoc assumption of current 
synthesis codes ( Shore! 20031. 



Stellar evolution theory, stellar atmospheres theory and the stellar birth 
rate are the basic input ingredients of evolutionary synthesis codes, which are 
tools designed to solve Eq. |^1 From the point of view of synthesis codes, these 
input ingredients can be classified in two categories: input data and input pa- 
rameters. Stellar tracks and model atmospheres, which are often included as 
built-in libraries, are the input data. The stellar birth rate and a few other 
quantities, such as the lower and upper mass limits of the IMF, the age and the 
metallicity, are the input parameters, that is they are switches that tell the code 
how to combine the input data to produce a model population. Usually, the 



1 Here we adopt the definition of stellar evolutionary phase used in stellar evolutionary theory, i.e. 
a particular stage of stellar evolution characterized by well-defined structural features, usually 
a given mode of nuclear burning: examples of evolutionary phases are the main sequence and 
the red giant branch. We warn, however, that the expression 'stellar evolutionary phase' is 
sometimes used more loosely in population synthesis literature, indicating any of the arbitrarily 
defined classes of Eq. [5] 
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input data available to a given code are selected and controlled by the author of 
the code. 

The general structure of synthesis codes is summarized in Fig. ^ Input 
ingredients will be described in the next section. 



Evolutionary tracks 

Metallicities 
Mass loss rates 



Stellar birth rate 

Initial mass function Star formation history 
Slope star formation rate 

Lower and upper mass limits Instantaneous burst 

Constant star formation 



Rotation 

Angle distribution 

Rotational velocities distribution 




Atmosphere libraries 



Evolutionary synthesis 



Stellar population 
(HR diagram) 



Synthetic spectra 
Colors 

Relative contributions 
Emission line spectrin? 



Figure 1. General structure of evolutionary synthesis codes. 



2.1. Evolutionary tracks and isochrones 

The stellar models entering a synthesis code are made up of two different parts: 
evolutionary tracks and stellar atmospheres templates. The former describe 
the inner structure and evolution of stars, and the latter their outer, radiation- 
emitting layers. The computation of these two parts require completely dif- 
ferent physics and is carried out by different codes: hence, evolutionary tracks 
and model atmospheres or empirical libraries are provided separately and then 
coupled together in the synthesis code. 

Evolutionary tracks are so called because they carry information on the 
path (track) that a star follows on the HR diagram throughout its lifetime (evo- 
lutionary). This naming convention is somewhat reductive, since evolutionary 
tracks in fact include much more information than the one displayed in the HR 
diagram: in fact, each time-point of an evolutionary track also contains de- 
tails of the stellar structure, that is the dependence on radius of variables such 
as temperature, density, chemical composition, nuclear energy production rate, 
etc. Therefore, they are in fact models of stellar structure and evolution: here, 
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however, we will follow the customary use and refer to them as evolutionary (or 
stellar) tracks. At first order, an evolutionary track is identified by its initial 
mass and metallicity. 

The tracks determine strongly the scope and the limitations of a synthesis 
model, so the choice of which tracks to use is an important one both for code 
developers and for end users. In the following, we will define a set of stellar 
tracks as a batch of stellar tracks of different initial masses computed under 
homogenous assumptions. Sets of stellar tracks available in the literature differ 
from each other in many respects: the metallicity values available, the mass 
range covered, the treatment of mass loss, the inclusion of rotation, plus other 
lesser-order effects such as the treatment of convection and the opacity tables 
included. Stellar track makers are often specialized in a small niche of stellar 
evolution theory, e.g. low-mass or massive stars. For this reason, the choice of 
the stellar tracks to be used in a synthesis model depends a lot on the kind of 
model one wants to build. For example, a synthesis model of an old population 
requires accurate tracks of low-mass stars, while a model of a young population 
requires accurate tracks of massive stars. 

Ideally, stellar tracks with any desired value of the stellar parameters should 
be available. Unfortunately, since evolutionary tracks are costly to compute 
and store, available models have only a very limited subset of all the possible 
metallicity and initial mass values. Typically, a given author usually provides 
sets of tracks corresponding to less than ten metallicity values to cover a range 
of several orders of magnitude. A second issue related to evolutionary tracks is 
that a given set of tracks does not necessarily cover all the evolutionary phases. 
In such a case the synthesis model maker must choose between keeping a self- 
consistent set of (incomplete) tracks or mixing different sets of tracks (sometimes 
with different physical assumptions). Therefore, when using or computing a 
synthesis model, one should be aware of the specific features of the stellar models 
included in it. 

Once the set of evolutionary tracks has been defined, one has a set of 
functions that describe the evolution as a function of time t of stars with 
given initial mass m, ftra,ck(t', m )- The next step is transforming it into an 
isochrone. An isochrone is defined by the population of coeval stars of differ- 
ent masses at a given age. Isochrones can be represented in different plans: 
as example, an isochrone in the theoretical HR diagram is a parametric curve 
/i SO (?7i;£) = [Lboi(?w; t), T e Q (m; t)] (m being the initial mass). 

Isochrones are like snapshots of evolving populations. Conceptually, they 
are an intermediate step in the transition from information on individual stars 
to information on composite populations. Not all codes perform this step, as 
isochrones are available in the literature so they can directly be used in a code, 
bypassing the need for evolutionary tracks. But those codes that use tracks as 
input must interpolate between evolutionary tracks to obtain the isochrones at 
a given value of t. 



2.2. Atmosphere libraries 

The atmosphere is the outermost part of a star. The mass of a star's atmosphere 
is negligible when compared to the star's total mass, so it has no effect on the 
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stellar structure and evolution; but, by definition, its optical depth at most of 
frequencies is small, so it is the region where the emitted spectrum forms. 

A synthesis code needs an atmosphere library to transform the theoretical 
quantities obtained with the evolutionary tracks into observational ones. There 
are two possible alternatives: 

1. Grids of theoretical atmospheres, composed of atmosphere models. The 
main challenge of the models is that in the atmosphere the radiation field 
is strongly anisotropic, so the transfer equation must be solved explic- 
itly. Different approximations can be adopted to solve this problem in 
model atmospheres, and one must be aware of the implications of choosing 
among different models, as different authors make different assumptions 
when they face the extraordinarily complex problem of computing a stel- 
lar atmosphere. Model atmospheres are generally given as a function of 
metallicity, gravity (g) and effective temperature. 

2. Grids of empirical atmospheres, composed of observed stars. These stars 
must be calibrated in flux and their physical properties must be well de- 
termined. Since the observed stars are used as templates for the stellar 
classes that have a bolometric luminosity different from the one of the 
observed star, a normalization it is also necessary. 

Analogously to what happens with evolutionary tracks, the available models 
cover only a very limited subset of all the parameter space. For a more detailed 
paper about the atmospheres libraries used in synthesis codes, we refer to the 
contribution by Gustavo Bruzual in these proceedings. 

2.3. Stellar birth rate 

In order to obtain the weight Wi of each stellar class, we need a recipe to compute 
the number of stars populating each class at any given time. To this scope, 
we need to know which kind of stars occupy a certain position at the time 
considered, and how many of them there are. The first piece of information is 
given by evolutionary tracks combined with model atmospheres, which we have 
introduced in ^i.l.l and H'2.2.1 The second piece of information is the stellar birth- 
rate. In order to simplify the subject, we will assume in the following that the 
mass and the time dependence of the stellar birth rate can be separated, as most 
synthesis codes also assume. This is a very strong hypothesis, which is assumed 
ad hoc to simplify the handling of the equations, and which is not necessarily 
a realistic one: for example, it is trivially false if stellar winds from massive 
stars inhibit the formation of low mass stars, or in any oth er case in which the 
IMF k e eps memory of pre v ious star formation epis odes. See lTenorio-Tagle et al.l 
1 200fiMShore et aU Jl987); S hore fc Ferrinll (|1995T ) for more details on this topic. 



For simplicity, let us also assume that the SFH is described by a Dirac delta 
function (the star- formation episode is extremely concentrated in time). This 
scenario has been labeled in different ways as coeval star formation, instanta- 
neous burst or simple stellar population (SSP). This mathematical approxima- 
tion allows us to neglect the chemical evolution of stars in the cluster (since all 
the stars are formed with the same initial metallicity) and to define a zero-point 
in the timeline. 
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To compute the expected emission of a stellar population in the SSP hy- 
pothesis, it is only necessary to provide its initial mass composition, which is 
defined by the IMF. The concept of IMF is rooted in measurements of the rel- 
ative frequencies of stars with different masses, which are observed to be fairly 
constant across different stellar populations. The IMF is often approximated 
either by a power-law or by a sum of power-laws over different subranges: 

. . dN _ a 
m ™ = 3— oc m . (3) 
dm 

A typical valu e for the power-law exponent is the one by Salp eter, a = 2.35 
(ISalnetenll955T) . Alternative IM Fs have also been proposed (e.g. iMiller fc Scalol 
119791 : lRanalll987l : iKroup J200lh . The IMF in the low-mass tail is poorly known, 
due to incomplete detection. The IMF in the high-mass tail is poorly known, 
due to the intrinsically low mass counts. We refer to these proceedings for more 
extensive reviews on the subject. 

2.4. Types of synthesis codes 

With these ingredients it is now possible to compute the contribution of the 
different stellar classes: since each stellar class present at a given time can be 
defined by a range of initial stellar masses, its contribution, iV{, is simply given 
by: 

Wi= tpu.(m) dm, (4) 

J m up 

where m l ° w and are the lower and upper limits of the z-th mass bin that 
gives rise to the i-th stellar class (the specific way in which these limits are 
defined varies from code to code). The codes that compute the contributions of 
the different stellar classes in this way, based on mass bins, are usually called 
isochrone synthesis codes. 

There is at least another way to compute these contributions, which is based 
on the evolutionary times of the different phases. To understand this point, let's 
consider how an isochrone of a given age is populated. The post-MS portion of an 
isochrone is populated by all those stars who have already evolved out of the MS 
but that have not died yet. Since post-MS evolutionary times are much shorter 
than MS times, the difference in mass between the turn-off (TO) star (the star 
that is just about to leave the MS) and the most evolved star is comparatively 
small. It can be shown that the older the SSP, the smaller such difference: 
asymptotically, all the post-MS isochrone portion is populated by stars with the 
TO mass, and the isochrone tend to converge to the corresponding evolutionary 
track to the degree that, from the TO on, they merge completely. This is shown 
with an example in Fig. [21 which illustrates the point in a less extreme situation 
(i.e. one in which the merging point takes place at a mass point located beyond 
the TO). 

The consequence of the merging of isochrones and evolutionary tracks for 
late evolutionary stages is that the post-MS luminosity of a SSP population is 
proportional to the fuel available to the TO star. Therefore, the contribution 
Wi of a given post-MS stellar classe can by expressed by Eq. |SJ 
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Sch aller et all iTtOOl . The mass point where the isochrone is equivalent to 
the track, 39 M Q , is marked in the diagram. 



Wi= (pM{m) dm = ipM(m cq ) Srrii = ipM(m cq ) 

Jm up 



dm. 



cq 



dt 



St; 



(5) 



where m eq is the mass point in the isochrone where the track and the isochrone 
become equivalent, and 5ti is the duration of the evolutionary segment corre- 
sponding to the i-th stellar class. This way of computing the pos t-MS contribu- 
tion to luminosity is known as the fuel consumption theorem fFCT: lRenzini k, Buzzoml 
1986; B uzzon i 1989). We refer to these papers for further details. Synthesis codes 
based on Eq. [5] are usually called fuel consumption codes. 

In the practice, the distinction between isochrone synthesis and fuel con- 
sumption codes is not a sharp one, as most codes make use of both methods. 
For example, since the FCT makes explicit reference to post-MS phases, FCT 
codes compute the contribution of the MS as isochrone synthesis codes do. On 
the other hand, some isochrone codes make use of the FCT to include post- 
Asymptotic Giant Branch (post-AGB) evolution. 



3. Overview of the uncertainties 



We can now move on to discuss the main sources of uncertainty in synthesis 
codes. These can be grouped as follows: 
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1. The uncertainties in the evolutionary tracks used and in the transformation 
from tracks to isochrones. 

2. The uncertainties in the assumed stellar birth-rate. 

3. The uncertainties in the assumed atmosphere libraries and their assigna- 
tion to theoretical quantities. 

4. The incompleteness of input ingredients. 

5. The numerical methods used. 

These uncertainties will be discussed in the following. Throughout this 
discussion, it will be useful to keep as a reference the sketch of the general 
structure of synthesis codes shown in Fig. ^ 

3.1. Uncertainties related to the evolutionary tracks and isochrones 

The homology assumption in track interpolations Some track providers (but 
not all) tabulate their tracks following a sequence of equivalent evolutionary 
points. These are defined as points in stellar tracks of different initial masses 
that can be related by means of homology relations: i.e. the structures of two 
models star of different ini tial masses are homologous a t corresponding equiy - 
alent evolutionary points (|Kippenhahn fe Weigertl Il99d : iMowlavi et alJ LL998). 
For exampe, in the case of a completely radiative main sequence (MS) star with 
a Kramer's opacity law the homology relations are: 



where kq is the opacity, [i the mean molecular weight, eo the energy production 
rate and X the fraction of hydrogen in the stellar core. These relations, together 
with the Stefan-Boltzmann law, 



allow to perform interpolations among tracks in the MS, assuming that the 
mass is constant throughout the evolution, thus obtaining tracks that are not 
tabulated. These relations provide an interpolation scheme to obtain £bo\( m ) 
and T e ff(m) for any desired mass. 

In fact, homology relatio ns are a good approximation only for MS stars 
( Kippenhahn fc Weigerttf l99C0 : beyond the MS the homology assumption always 
breaks down. However, in isochrone computations the interpolation scheme 
implicitly assumes homology even beyond the MS. 

A second point to be addressed is that due to the presence of mass loss the 
assumed homology relations are not valid, since mass loss is coupled with the 
stellar evolution in a non-trivial way that is not considered in the usual homology 
relations; but, again, the interpolation scheme among tracks is maintained even 
if tracks with mass loss are used. 




(6) 




(7) 
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Summing up, the interpolation scheme which is customarily used to cal- 
culate tracks that are not tabulated is based on an assumption which is, most 
of the times, not valid. There are, however, some ages at which the homology 
assumption has only a minor impact on the computed models. These are the 
ages at which the stars with tabulated tracks are at the TO. As we have shown 
before, post-MS stars have masses similar to the TO mass; if the TO mass is 
a tabulated one, the contribution of post-MS stars will not be much affected 
by interpolation errors. This observation permits to formulate a criterion for a 
reliable comparison among synthesis models: 

The most reliable ages for model use or comparison are the MS turn-off ages of 

the input tracks. 

If differences appear at such ages, they reflect differences in the numerical 
methods used in each code, which must be further investigated. This criterion 
is only valid, of course, if the codes have the same input ingredients. 

Fast evolutionary phases In post MS phases, luminosity is not always a well- 
behaved function of the mass, so the isochrone is not always defined: a typical 
isochrone features shallow sections as well as peaks and discontinuities in the 
m,£ plane, as shown in Fig. |31 Shallow sections correspond to quiescent phases 
of stellar evolution, where evolution is slow (e.g. the MS); peaks correspond to 
faster phases (e.g. the asymptotic giant branch); and discontinuities correspond 
to abrupt transitions between phases (e.g. the onset of central helium burn- 
ing in intermediate mass stars) or jumps in stellar behavior (e.g., the transition 
between Wolf-Rayet [WR] and non-WR-type structures). Peaks and discontinu- 
ities will generically be labeled in the following as fast evolutionary phases. It is 
interesting to realize that the expression fast evolutionary phases makes explicit 
reference to time variation, although by definition time is a constant throughout 
an isochrone! In fact, the term fast refers to the behaviour of the star whose evo- 
lutionary track coincides with the isochrone in the post-MS phases. In terms of 
isochrones, the expression refers to regions where the derivative di/dm diverges. 

Fast evolutionary phases are difficult to handle, because a small difference 
in the initial stellar mass yields a large difference (or an indetermination) in 
luminosity, so that the result depends crucially on which luminosity is chosen 
to be representative of a given stellar class. In principle, this difficulty can be 
dealt with by defining stellar classes so that they are characterized by small lu- 
minosity ranges and they do not go across discontinuities; however, the available 
resolution in mass is governed by the evolutionary tracks used by the code, and 
is generally too low to resolve adequately such phases in synthesis models. This 
is a severe problem, since fast evolutionary phases are ubiquitous in post MS 
evolution, and at certain frequencies they bear a major weight in the luminosity 
balance. 

This problem has been present in synthesis codes since its very beginning 
( Tinslev Sz Gunnlll976f ). In the particular case of discontinuities, the way in 
which it is tackled is often labeled a 'technical detail' of the computation and 
dismissed as unimportant, and thus the papers describing evolutionary synthesis 
models do not generally make any reference to its solution - in spite of its dif- 
ficulty and of the potentially disastrous consequences of incorrect assumptions. 




Figure 3. Details of fast evolutionary phases in the V (solid line, left axis) 
and K (dotted-line, right axis) bands. Bottom panel: complete isochrones. 
Middle panels: blow-up of the mass axis, showing details of fast evolutionary 
phases at three different ages. Top panels: same as middle panels, with a 
more extreme blow-up. The isochrones have been computed bv l Girardi et all 
(2002J from the evolutionary tracks bv lMarigo fc Girardil <|2001f) 



Here are a few examples of the ways in which the pro blem of discontinuities has 
been approached: a) in the Starburst99 synthesis code (|Leitherer et al . 1999), for 
certain metallicity values, an undocu mented stellar t r ack a t 1.70 1 M m is added 
to the tabulated track of 1.70 M Q bv ISchaller et al.l ( 19921 ) and lSchaerer et all 
( 1993al ). to avoid stellar classes going across the discontinuity of the stellar 
models' behavior at such mass (C. Leitherer, D. Schaerer, & G. Meynet, pri- 
vate communication) ; b) to deal with the same problem, additional evolutionary 
tr acks around the same m ass ra nge are used in the co mputation of the isochrones 
bv ICastellani et aD (I2003T) and ICariulo et al.l (|2004l ) (S. Degl'Innocenti, private 
communication) and lBrocato et alJ (|l999h (E. Brocato, private communication); 
c) to avoid the intrinsic di scontinuity in t he iso chrones, the same mass is used 
twice in the isochrones by iGirardi et al.l (|2002l ) , namely at the end of the red 
giant branch (RGB) and at the beginning of the horizontal branch (S. Bressan, 
private communication) . 

This problem shows up in the results of synthesis codes. For example, Fig. 
0] in which the number of post-MS stars is plotted as a function of time, shows a 
gap at ages around 10 9 years; the gap corresponds to the MS evolutionary time 
of stars with initial masses near the 1.7 M discontinuity. 
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Figure 4. Number of post MS stars as a functio n of time for different metal- 
licities, obtained from the isochrones compiled bv lGirardi et alJ f|2002f) . 



The figure also shows the presence of small peaks at young ages. The 
maximum of each peak indeed corresponds to the MS lifetime of the tabulated 
tracks. The peaks appear as a consequence of an incorrect time interpolation 
scheme and can be corrected by adopting a different one. This problem will be 
discussed in the next paragraph. 



Evolutionary times and time interpolations A further problem related to the 
computation of isochrones is the duration of the different evolutionary phases 
that are included in the tracks. In ^12.1.l we have shown that, in order to obtain 
isochrones, it is necessary to interpolate among the tabulated tracks so that 
tracks with intermediate values are obtained. After that, it is necessary to in- 
terpolate in time among different tracks to obtain the isochrone at a given age. 
The usual scheme assumes lineal interpolations in a logarithmic space since, 
at first order, the duration of an evolutionary phase is proportional to the ratio 
between mass (available fuel) and luminosity (consumption rate). Under the ho- 
mology assumption, these are related by a power-law dependence. Again, such 
an approximation may not be a good one if homology fails. The potential con- 
sequences are clearly illustrated by Fig. 03 which shows the supernova (SN) rate 
obtained using different interpolation schemes. The saw-teeth behavior of the 
solid line is clearly an arti fact of the line a r inte rpolation scheme. The parabolic 
interpolation proposed bv I Cervino et alJ (|200lT ) yields a much smoother behav- 
ior. 



Uncertainties in population synthesis 



13 



^3 
CO 



CD 

cd 



2 CO 



<IOMo 



25Mo 20Mo 



15Uo 



12Mo 



85Mt> 



l20Mo 



v 



60Vc 




Lnco r 

Parabolic 



6.4 



7.2 



Figure 5. SN rate vs. time using different interpolation techniques. The 
solid line corresponds to a linear interpolation in the logm — logt plane. The 
short-dashed line corresponds to a parabolic interpolati on. The evolutionary 
tracks used have been computed bv lSchaller et alJ 1)199: 
insta ntaneous; solid line, a = 2.35, M up = 120 M Q ; i 
from lCervmo et alJ (|2001fl . 



Star formation law: 
0.020. Figure taken 



The accuracy of interpolation schemes for stellar evolutionary times has 
barely been address ed in the literature. In fact, the solution proposed by 
ICervino et alJ ( 200lh is only a second order mathematical approximation, that 



does not take the physics of stellar evolution into account. A more accurate 
estimation of the lifetime of different evolutionary phases can be obtained by 
fitting the evolutionary times of different evolutionary phases. As an example, 
we show in Fig H3 the MS lifetime as a function of the initial mass for different 
sets of evoluti onary tracks. The fit of all the points of various sets of tracks 
( Buzzonill2002t l provides a dependence of the MS lifetime with mass: 



log t = 0.825 log 2 + 6.43. (8) 
120 w 

This simple formula shows how the interp olation scheme of parabolic interpola- 
tions can be greatly improved. We refer to Buzzoni (2002) for additional details. 

Note that any variation in the time interpolation scheme not only affects 
the SN rate, but also any other predicted quantity that depend on lifetimes. In 
particular, the stellar make-up of the synthetic cluster also changes. 
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Figure 6. Theoretical MS lifetime vs. stel lar mas s for s o lar metallic- 
itv ac c ording to the e volutio nary tracks bv iBeckerl ll 19811): Vandenberel 
JlQS.qft: ICastellani et all i|1990D : iLattanziol l|1991[k ISchaller et all l|1992ft and 
iBressan et alJ (11 9931). The solid line is a fit to the data according to Eq. (3). 
Figure from Buzzoni (2002). 



Further inconsistencies in the computed isochrones As we have previously seen, 
all types of synthesis codes rely on interpolations among tabulated evolutionary 
tracks independently of the method they use. Usually, these interpolations are 
linear in a logarithmic space for quantities which are assumed to follow homology 
relations, such as luminosity or effective temperature. This interpolation scheme 
is applied, without any precise physical reason, to other relevant quantities also: 
for example, abundance ratios in stellar atmospheres (needed to determine the 
presence of WR stars in young populations), or the mass loss rates (needed to 
compute the amount of kinetic energy deposed in the interstellar medium). 

To illustrate the spurious effects of such ad hoc interpolation schemes, we 
show in Fig. [7|the ratio R between the mass ejected by a star as obtained from 
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Figure 7. Ratio of the integrated mass-loss during the lifetime of the star 
vs. "reconstructed" mass loss by subtraction o f the mass at the the end 
of the evolution for the tracks bv ISchaller et alJ l|1992l) (solid line) at solar 
metallicity and standard mass- loss rate and the tracks bv lMevnet et al.l |l994) 
at twice solar metal licity and twice mass- loss rate tracks (dashed line). See 
iCerviho et a D il200l for more details. 



integration of the mass loss rate, and the difference be tween the initial and final 
masses as given by the track (see ICervino et alJl200li for more details): 

R = Eti"fa)fe-^-i) (Q) 

M-m ke { } 

where the index k e refers to the last tabulated point in the track. Consistency 
would require such a ratio to be equal to 1. However, R is found to take values 
of up to 1.4, i.e. the integrated mass lost by stars is inconsistent with the 
'structural' value by up to 40%. 

3.2. Uncertainties related to atmosphere libraries 

As mentioned in ^2.2.1 stellar libraries are needed to transform the theoretical 
space (T e g, g and () into the observational one (colors and/or spectral energy 
distributions). Usually, the available model atmospheres form a coarse grid in 
the (log g, log T e ff) plane, whereas isochrones are generally continuous in the 
plane. In order to assign a model atmosphere to each isochrone location, one can 
either choose the nearest atmosphere of the grid (closest-model approximation), 
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or interpolate between nearby atmospheres. Assigning the nearest atmosphere 
implies a further binning of data and may originate jumps in the results. 

54 I — I — I — I — I 1 — I — I — I — I 1 — I — I — I — I 1 
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log (Time [yr]) 

Figure 8. Number of ionizing photons as a function of age for solar metallic- 
ity and different IMF slopes. Note the downstairs-like behavior of the quantity 
at evolved ages due to the transition among atmosphere models. Figure from 
iLeitherer etafl l|1999fl . 

It is often assumed that these jumps cancel on average when a whole popu- 
lation is considered. Unfortunately, this is not always the case, as is made clear 
by Fig. EJ which shows the evolution of the number of ionizing photons in a 
cluster as a function of age. The downstairs-like behavior of the plotted quan- 
tity at evolved ages is due to the transition among atmosphere models arising 
from the use of the closest-model approximation. 

An alternative solution is interpolating among atmosphere models to ob- 
tain smoother results. This solution is not common in synthesis codes, since 
it increases the computational time. However, it is mandatory in the analysis 
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of color-magnitude diagrams (CMD), and, as we will see later, it will also be 
mandatory in synthesis models including rotation. Therefore, we mention here 
some of the unsolved issues related to this kind of interpolations. 

If the interpolation scheme is based on the Stefan-Boltzmann law (Eq. 0) 
the interpolations are wavelength independent, and it is simple to perform them 
linearly in the log^ — logT c g- plane. However, for the case of a spectral en- 
ergy distribution it is more realistic to consider the emission as a black-body 
spectrum: 



logoff) = logf^) -h»o(r' - J 



V A 5 

= tog^^^-Or + logCl-e-*), (10) 

where x = hc/XkT e g. This means that, at a given wavelength, log^\ is a linear 
function of the variable y = log(e 1 ' — 1) = x + log(l — e~ x ). Such interpolation 
scheme clearly differs from the one based on the Stefan-Boltzmann law. In 
the intermediate case, i.e. the computations of broad or narrow bands, the 
interpolation scheme should be flexible enough to tend towards both extreme 
schemes (Stefan-Boltzmann law and black- body law) de p end ing on the width 
of the band. We refer to the appendix in iJamet et al. for a detailed 

description of this issue. 

3.3. Incompleteness of the input ingredients 

Incompleteness of evolutionary tracks The lack of homogeneous sets of stel- 
lar models that cover all the evolutionary phases is perhaps the most severe 
problem in population synthesis. When a synthesis model is to be applied to 
observations, the lack of homogeneous stellar models must be sidestepped in 
some way, sometimes sacrificing homogeneity, sometimes using inputs based on 
physics that is known to be incorrect. Here are some examples: (1) The inclu- 
sion of the post-AGB phase in synthesis codes requires using tracks computed by 
different authors (with possibly different input physics), and performing some 
ad hoc assumptions to link the post-AGB to the AGB tracks. (2) Stellar models 
computed with enhanced mass loss rates are assumed ad hoc to mimic the ef- 
fects of rotation in massive stellar evolution, although observations point towards 
mass- loss rates even lower than the standard ones. 

In such cases, assessing the uncertainty in the results is next to impossible, 
particularly due to the lack of a physical picture giving a realistic guide of 
the problem. However, these are the only way of producing a complete result 
when the models are computed for comparison with the observations, since they 
include all the relevant phases of stellar evolution. On the other hand, when the 
goal of a synthesis model is to link our (incomplete) knowledge on stellar theory 
to the emission from an ensemble of stars, incomplete input data can be used. 
In this case, the results will be incomplete, but they will at least predict a lower 
limit to the expected emission. 



Incompleteness of atmosphere templates The region covered by the isochrones 
in the log g — log T e g plane is not the same as the one covered by the atmosphere 
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Figure 9. Right: log q vs. logT e gf coverage of the ST ELIB empirical li- 
brary l|Le Borgne et al.ll2003ft used bv lBruzual fc Charlotl (2003). The differ- 
ent symbols represent different stellar spectral classes: full circles are dwarf 
MS stars (class V), open circles are giants (class III), and plus signs are su- 
pergiants of classes I and II. Small circles are used for stars with no spectral 
class determination. Left: Stellar library (large symbols) in the plane logg vs. 
logToff, comp ared to the position of stellar models of solar me tallicity (small 
dots) used bv lGirar di et alJ II2002D. The spect ra are taken f r omlCastelli et"aTI 
<|1997|) ( crosses') . iFluks et alJ (|1994f) (circles') lAllard et all <|2000j) (squares'! . 
and pure blackbody (tria ngles). (See lGirardi et al.ll2002l for details). Figure 
from lCdrardi"et~aTl J200l . 



libraries. Specifically, theoretical libraries do not cover some important regions of 
this plane, like the one occupied by stars with T e g around 10 4 K and intermediate 
log g values, that corresponds to the A supergiant region, and the one occupied 
by cool stars (T c g < 3000 K) with negative log g values, that corresponds to the 
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red supergiant region. To bypass this problem, some codes make use of empirical 
libraries. 

A severe limitation of empirical libraries is that they do not have a good 
coverage in metallicity, a fact which limits the possibility of applying the models. 
Furthermore, although the red supergiant region is well populated (since these 
are luminous stars) , the coverage in other regions of the log g — log T c g plane is 
not sufficient to produce realistic results even at solar metallicity. 

In Fig. El the different areas of the log g vs. log T e g- diagram covered by 
empirical and theoretical libraries can be seen: note how neither of them covers 
the whole diagram. Given these limitations, it is impossible for a code to produce 
a physical self-consistent result. 

Inadequate sampling of the evolutionary phases A further problem for the com- 
putation of synthesis models is the insufficient sampling of the stellar evolution 
along the tracks. This problem has a different impact depending on the scope 
of the synthesis model and the way it is compared to the observational data. In 
the case of population synthesis, one wants to reproduce the total luminosity 
of the cluster based on the integrated emission, so one needs to take into ac- 
count carefully the emission of all the stars. On the other hand, when studying 
CMDs one wants to reproduce the emission of the stellar population based on 
its representation in the HR diagram; hence, not all phases are necessary. This 
difference has the following consequences: 

• In the case of CMD it is not necessary to cover all the HR diagram, but 
only the areas under study. Relevant points in such areas should be in- 
cluded explicitly in the isochrones since they can be directly compared 
with individual stars. This is the case of the TO, the tip of the RGB or 
the tip of the AGB. These points are usually included in the evolutionary 
tracks and the isochrones. 

• In the case of population synthesis the objective is to obtain the integrated 
light from the ensemble. So suitable isochrones should be defined by points 
that are a correct representation of all the evolutionary phases, which 
are not necessarily the extreme points. For example, if a stellar class is 
defined around the mass point of maximum luminosity, the corresponding 
contribution is overestimated. 

• Not all the stellar models computed to follow the stellar evolution are 
eventually published (Georges Meynet and Daniel Schaerer, private com- 
munication; Antonio Claret, private communication). The reason is that, 
in regions where real stars are barely observed, the advantage of includ- 
ing these points is not considered worth the extra-memory needed. Hence, 
evolutionary tracks are reduced to a minimal set of equivalent evolutionary 
points plus a few additional points that are relevant for stellar evolution 
theory. However, population synthesis would require the detailed knowl- 
edge of the evolutionary path in the HR diagram. 

Summing up, evolutionary tracks are optimized for the comparison with 
CMD (which provides a first reliable check of the assumed physics), but not to 
perform evolutionary synthesis studies. These difference is clear when extreme 
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phases of stellar evolution are considered. For example, the tip of the RGB is a 
useful point for representation in CMD, but it is not a suitable point for pop- 
ulation synthesis, since it represents a transient high luminous phase in stellar 
evolution that pollutes the emission of the ensemble. As a further example, SN 
explosions are not included in synthesis models because they dominate the emis- 
sion with their high luminosity (i.e. the resulting model would have a spectral 
energy distribution that corresponds to a single component, the SN explosion, 
whatever the underlying population is). 

The previous conclusion can be restated in a different way: since synthesis 
models are designed to obtain the emission of an ensemble of stars, and not 
individual stars, their use should be limited to situations where an ensemble 
rea lly exists. A necessary condi tion to fulfil this requirement has been formulated 
by ICerviho &: Luridianal ( 2004T ) by defining the Lowest Luminosity Limit that 



an observed cluster must have in order to be modelizable with synthesis codes: 

The total luminosity of the observed cluster must be larger than the individual 
contribution of any of the stars included in the model. 

Quite obviously, the inclusion of SN events in the model increases the cor- 
responding Lowest Luminosity Limit. 



4. Uncertainties related to the numerical methods 

The final source of uncertainty in synthesis models lies in the numerical methods 
used. Three types of uncertainties arise at this level: the presence of bugs in 
the program, the accuracy of the numerical methods, and the consistency of the 
algorithm used. 

4.1. Bugs 



As observed by Ferland (http://www.nublado.org), the presence of bugs is in- 



evitable in any large code. Sometimes they are documented and solved, and 
corrected in following versions. As an example, we show here the SN rate ob- 
tained by the original version of Starburst99 (c.f. Fig HOj) . The figure can be 
compared to Fig. |SJ which is a correct representation of the SN rate. Note in 
particular that the wrong SN rate increases with time, whereas the corrected 
one tends to decrease. 

The bug was fixed in Starburst99 following the prescriptions bv lCerviho et all 

The source of this bug was an incorrect determination of the m- ow and 



m^ p limits in Eq. ^ Although the bug showed up in the SN rate, it also affected 
the Wi determinations of all the other evolutionary phases. The effect of the 
bug was particularly evident in the case of the SN rate because, being this a 
derivative, it amplifies any discontinuity in the primitive function. 

4.2. Accuracy 

The second source of uncertainty is the accuracy of the numerical methods used. 
As an exampl e we p lot in Fig. ^2 the SN rate of Type I SN computed by 
iRomano etahl (|200fih . In this case, the plot shows a correct overall behavior. 
However, at age larger than 10 Ga, an oscillation of the results around a mean 
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Figure 10. SN rate vs. time obtained by the original version of Starburst99 
code. Star formation law: instantaneous; solid line, a = 2.35, M up = 100 M ; 
long-dashed line, a — 3.30, M up = 100 M Q ; short-dashed line, a — 2.35, M up 
= 30 M Q ; Z = 0.020. This figure can be compared with Fig. Efl The bug was 
fixed i n Starburst99 code following the prescriptions shown in lCerviho et al.l 



value is observed. This behavior is caused by insufficient numerical precision. 
This oscillation had only a minor effect on the original paper, which deals with 
chemical evolution, a field in which the relevant quantities come from cumulative 
contributions; but if this model were applied to stellar population synthesis, the 
nominal tabulated result could not be used. 
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Figure 11. Ty pe la SN rates obtained w ith different assu mptions f or the 
stellar lifetimes: iMaeder fc Mevneij (119891) (thin solid line) : iTinslevI £l980) 
(thick dotted line): iSchaller et all ljl992j) (thick solid line); iKodamal l|l997|) 
(thick dashed line). The effect of changing the fraction of mass entering the 
formation of type la SN progenitors, A, is also shown. Th e type la SN rate 
observed i n the Galaxy a t the p resent time is also shown llCappell aro et ah 
199^ See|R omano etall lE)0l for more details. Figure from iRomano" et"ahl 



4.3. Algorithms 

FCT and isochrone synthesis codes have bee n shown to produce systematically 
discrepant results I Chariot &: Bruzual Il99ll ). As explained in ^2.4.1 the two 
types of codes use different algorithms to compute the contributions of the dif- 
ferent stellar classes. However, it must be emphasized that both types of codes 
are based on the same underlying physics. Furthermore, t he FCT can be taken 
into account in the isochrone com putations, as shown hv l Rressan et al. l jTijjfl 
and, especially. iMarigo k, Girardl ((200 ll ). Hence, both methods should produce 
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exactly the same results, provided they have, of course, the same inputs. Un- 
fortunately this is not the general case, a fact that suggests that the differences 
among synthesis codes computations could be due to differences in the interpo- 
lations scheme. Note that the current market of synthesis codes is dominated by 
isochrone synthesis codes that compare their results with similar models. Such 
comparisons do not shed any light on this problem, but only provide a con- 
sistency test for the mathematical methods used by isochrone synthesis. Only 
comparisons among different methods will increase the reliability of evolutionary 
synthesis: the stellar population synthesis method will not be a reliable tool until 
the isochrone and the FCT methods yield consistent results. Although this may 
seem a pessimistic point of view, especially to those who make use of synthesis 
codes to infer the physical properties of observed clusters, the situation is not so 
dramatic: in spite of persisting differences in the results, there are also regions of 
the electroma gnetic spectrum where both methods provide similar results (see 
lBuzzonJl2005l . and Fig. ^Jfor a comparison). Hence, such wavelength regions 
can be confidently compared to observational data. 

4.4. Documentation, documentation, documentation 

In this section we have shown that, besides the possible uncertainties in the in- 
put ingredients, there are also uncertainties associated to the methods used. The 
evaluation of such uncertainties is quite difficult since they depend on the code 
itself. The only feasible way to evaluate how reliable the model results are is to 
write an exhaustive documentation of the corresponding code, both at the level 
of a user's guide and at a technical level. A good example to follow is the docu- 
mentation of the photoionization code CLOUDY (http://www.nublado.org). 

5. Rotation and variability 

In the previous sections we have addressed the uncertainties in synthesis models 
in the simplest and most idealized case which is the assumption of the existence 
of a single isochrone that defines the stellar mixture at a given age. But in a 
stellar cluster it is possible to have stars with distributed features, such as stars 
that rotate with different rotational velocities, variable stars or binaries, such 
that their overall emission is also distributed. In terms of modeling, this means 
that a well-defined spectral energy distribution (the main output of synthesis 
models) is not enough to characterize the cluster, and that a distribution range 
must also be provided. 

The case of rotation is illustrated in Fig. I13l where a set of isochrones at the 
same age is plotted. Rotation produces a loss of symmetry in the star, so that its 
emission depends on the inclination angle of the star with respect to the observer. 
Points in Fig. ^] differ from each other in the initia l mass, the rotational veloc- 
ity an d the inclination angle (see IClaretll2000l . 120031 : IClaret h Perez Hernande"3 
2005, for more details). 

The inclusion of rotating stars in synthesis models will imply a global revi- 
sion of synthesis models that includes: 

• Determining whether homology relations exist for rotating stars, in order 
to interpolate tracks and compute isochrones. 
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Figure 12 . Upper plot s of each panel: luminosity evolution of the SSP 
models bv lBuzzonil l|200fij) (solid dots), for solar metalli citv and Salpeter IMF 
compared with other theoretical outputs accordin g tolLeitherer et a~ ( 19991) 
(open dots') iBressan et all l|1994|) (dashed line) and B ruzual fc Chariot! |2003) 
(solid line). The total mass is scaled to M,c;cip= 10 1 1 Mm throughout, with 
stars in the range 0.1-120 M Q . The iLeitherer et all l)1999ft model has been 
slightly increased in luminosity (by about 0.06 mag in bolometric luminosity, 
at 1 Ga) to account for the luminosity contribution of low-MS stars (M* < 1.0 
M Q ). Lower plots of e ach panel : mod el residuals with respect to SSP fitting- 
functions. Figure from iBuzzonil <)2005|) . 



• Studying the distribution of initial rotational velocities and including it 
in the stellar birth rate. This also means making additional assumptions 
about whether the rotational velocity distribution and the IMF are sepa- 
rated distribution. 

• Including an inclination angle distribution in the stellar birth rate. Such 
a distribution should be a random (flat) distribution. 

• Abandoning the closest atmosphere approach in favor of an atmosphere 
interpolation scheme, and studying the variation of line profiles with ro- 



Uncertainties in population synthesis 



25 




(b-y) 



Figure 13. HR diagram for isochrones at the same age with different rota- 
tional velocities and inclinations angles. Courtesy of Antonio Claret 



tation (for high-resolution atmosphere templates): all the advantages of 
including rotation in the isochrones will be lost if, at the end, the combina- 
tions of spectral energy distributions do not take into account the richness 
of available situations and the isochrone points are described by a discrete 
set of atmosphere templates. 

• Obtaining not only average results but also the uncertainty associated to 
the mixture of distribution angles in the cluster. 



6. What is really computed? 

Although all evolutionary synthesis codes pursue essentially the same goal, they 
may differ substantially from each other in several aspects. Here, we will discuss 
two major distinctions within the field, the first related to the specific procedure 
followed in populating the IMF and the second to the interpretation of results. 

6.1. Modeling philosophies 

Evolutionary synthesis models can be either probabilistic or deterministic in two 
distinct aspects: a) the IMF sampling, and b) the interpretation of results. 



26 



M. Cervino and V. Luridiana 



The IMF sampling In a synthesis code, the IMF can be used either as an exact 
or as a probabilistic description of the mass distribution. In the first case, the 
relative frequencies of newborn stars of different masses are univocally fixed by 
the IMF; in the second, the mass of each newborn star is assigned through a 
random selection process where the IMF is the weighting function. In the first 
case, if the IMF is continuous, a smooth (although binned) distribution of stars 
results; in the second, the resulting distribution will be noisy, and the smaller 
the number of stars in the stellar population, the noisier the distribution. Two 
runs of the code with the same input parameters will yield two identical models 
in the first case, but not in the second (Fig. I14|) . The first approach is adopted 
by the so-called standard codes while the second is adopted by Monte Carlo 
codes. 

Standard codes are sometimes called analytical, but this adjective is mis- 
leading, because: a) in either case the underlying IMF may be an analytical 
function, and b) although the IMF is considered analytical, a discretization step 
(binning) is required by both methods. Hence, both methods are analytical in 
the same measure, and we will use here the less ambiguous expression standard. 

The interpretation of results As mentioned before, the main final product of 
a synthesis model is the expected value of the total luminosity. What does 
'expected' mean? Let's see here which are the possible answers. 

The interpretation of the output of Monte Carlo codes is straightforward: 
the luminosity computed by the code is the particular luminosity of a particu- 
lar realization of all the possible clusters characterized by a given set of input 
parameters. Running a set of Monte Carlo models with the same set of input 
parameters and data will produce a distribution of output luminosities. This 
distribution is a sampling distribution of the underlying Integrated Luminosity 
Distribution Function. The more numerous the stars in each simulation, the 
narrower the distribution (in relative terms). An example of this trend can be 
seen in Fig. 1151 The solid line in the middle of the colored areas is the prediction 
of a standard model at different wavelengths; the colored bands are the areas 
spanned by sets of Monte Carlo models with the indicated number of stars. A 
striking feature of this plot is that the uncertainty area depends critically on the 
wavelength considered. 

Another example is found in Fig. 1161 where the U— B and B— V colors re- 
sulting from different simulations in which stochastic fluctuations in the number 
of stars that populate different evolutionary stages is plotted with small dots. 
Lines corresponds to analytical simulations and me dium dots to LMC c l usters . 
The models have been computed with the code by iBruzual fc Charlotl (2003). 
From this figure it can be seen that Monte Carlo simulations for small clusters 
lie in a region that is not covered by the analytical simulations. This bias in 
colors is a natural effect when the integrated light coming from an observ ation 
with a low number of stars is analyzed (see ICerviho' fc Valls-Cabaudl l2003. for 
more details). 

In the case of standard models, however, different meanings can be attached 
to the output luminosity. Each time we run the code with a given set of input 
parameters we obtain the same result, so in this sense standard models are de- 
terministic. Indeed, the traditional interpretation of standard synthesis models 
is wholly deterministic: the output is interpreted as if it were the luminosity of 
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120 

Figure 14. IMFs computed with a slope a = 2.35 for 10 3 stars in the mass 
range 2-120 M©. It is compared here three stochasti cal Monte Carlo simu- 
lation s with an analytical function. Figure taken form lCervino fc Mas-Hessd 
llT99l . 

an (ideal) cluster with the same characteristic of the model. Models of this type, 
with their associated interpretation, are called deterministic synthesis models. 

Such interpretation is based on the assumption that the IMF is well sam- 
pled, which in turn implies that the IMF is seen as an exact description of the 
stellar population. This assumption is questionable, since whatever the mech- 
anism of star formation is, discrete entities cannot map continuously a curve. 
Even if we divide the mass range in bins, a continuous distribution would even- 
tually imply fractional numbers of stars in certain bins, which is unphysical. On 
the other hand, if we allow for a certain degree of stochasticity in the process of 
star formation (which, on physical grounds, seems plausible), these contradic- 
tions are overcome. But allowing for stochasticity in the IMF implies that the 
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Figure 15. Comparison of predicted luminosity of a standard model at dif- 
ferent wavelengths with the 90% confidence interval predicted by Monte Carlo 
models of 2. 10 5 M© (blue), 3. 1 4 M ffl (green) and 3. 10 3 M_odot (red). 
Models compu ted with sed@ code llCervino et alJl2000t ICerviho et alJ l2002b: 
ICervifio et al] f20021. Figure from ICervihol 1)2003') 



overall cluster properties will be stochastically distributed. This view is the one 
held by statistical synthesis codes. 

In sum, deterministic synthesis codes interpret the (deterministically ob- 
tained) output luminosity as the total luminosity of the cluster, while statistical 
synthesis codes interpret it as the mean value of a distribution of possible lu- 
minosities. The distinction between the concept of total luminosity and the 
concept of mean value of a luminosity distribution is a fundamental one: as an 
example, the mean of a distribution is not necessarily a value that the variable 
can take. Furthermore, knowledge of the mean value is not enough to estimate 
even roughly the variable's value, if the shape of the distribution is not known: 
that is, the mean value is neither an actual value of the total luminosity of a 
cluster, nor necessarily a good guess of it. 

Statistical and deterministic models are not actually different classes of 
models, but rather different interpretations of the same class of models, the 
standard ones. Some standard codes do not explore this interpretation and thus 
produce purely deterministic models, while others have built-in statistical tools 
that add to the main result, the luminosity, estimates of the other statistical 
parameters of the distribution. 

6.2. Statistical approximations to the problem 

The issue of sampling and the interpretation of the results of synthesis models 
has been addressed in the literature from different points of view. They include: 



Monte Carlo simulations and their comparison to standards models: these 
works make a description of the problem, but they do not provide any gen- 
eral solution other than performing tailored Monte Carlo simulat i ons fo r 
the cluster under c o nsideration. Examples a re: iBarbaro k, Bertellil (Il977f): 
IChiosi et ail(ll98Sf ): ISantos fc Frogell (|l997f t:lc erviho. Luridiana. Sz Castanderl 
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Figure 16. Comparison of Monte Carlo simulations (small dots), analytical 
simulations (lines) and observations of LMC clus ters (medium dots). Mod- 
els h as been performed b y the code presented bv lBruzu al fc C hariot! 1 2003) 
with iFaeotto et alJ l|1994j) evol utionary model s and lLeieune et al.1 l)1998j) at- 
mosphere models. Figure from iBruzuall (|2002f) . 



2nOCh lLancon fc Mouhcine] 12000): BruzuaJ (|2n02l 'l: lGirardil $2000): C antiello et all 



2003) 



Theoretical approaches based on the IMF sampling and Surface Brightness 
Fluctuations studies: These works try to quantify the impact of sampling 
effects on the interpretation of the results. Although relevant as a first 
step towards a characterization of the underlying luminosity functions, 
they are of limited power when it comes to obtaining the variance of the 
corresponding distribution. In this specific subject there have been two 
different approaches: 
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1. The variance of the luminosity function has been computed in terms of 
an effective number of stars TV, defined as the ratio between t he square 
of th e mean and the variance of the luminosity function ((Buzzonil 
Il989l ): 

1 _a 2 _ E Wjlf 

Af f (E Wik) 2 ' 1 ' 

This expression is used in the works bvlBuzzonil rtl989h:ICerviho et alJ 

(|2002bD : ICerviho fc Valls-Gabaudl (|2003l ): lGonzalez et al.M2004MCerviho fc Luridianal 
(|2004l ). and their use is focused on the determination of the uncer- 
tainty of synthesis models. 

2. The variance is presented in terms of Surface Brightness Fluctuations, 
L, defined as the r atio between the variance and the mean of the 
luminosity function ( Tonrv Sz Schneiderlll988l ): 



A* E Wik 



SBF are used to obtain extragalactic distances, and they have also 
been proposed as a means to break the age-metal licity deg e nerac y 
in old stellar p o pulations. They are used, e . g., by IWorthevl (Il994f); 
iBuzzonil (H993Tl: iLiu et all (l2000l ) ; iMei et al.l ipOOlfl : iBlakeslee et alJ 



2001); iGonzalez et alJ (120041 ): see also the contributions by Rosa A 
Gonzalez and Maurizio Salaris in these proceedings. 



AsH uzzonil ((19931 ) shows, both quantities are different ways to express the 



same quantity, and they express the statistical entropy of synthesis models 
(see also this proceedings or astro-ph/0509602 ). Given that both Af and I 
are estimates of the mean and the variance of the luminosity function, these 
quantities provides a powerful tool for the comparison of s ynthesis models . 
The differences in the I value found by different authors fsee lLiu et al . 2000; 
IMei et aDl200ll . for a comparison) suggest either a difference among the 
luminosity functions implicitly computed by synthesis models or problems 
in the numerical computations. 

Theoretical approaches based on luminosity functions. In other works, the 
luminosity function or its higher-order moments are studied. The approach 
is the o nly one that can prov i de confidence intervals of results. Exam- 
ples a re iGilfanov et al.1 ( 20041 ) ; ICerviho. Luridiana. &, Cerviho-Luridianal 



( 20051 ). 



Although the statistical modeling of stellar populations is not extensively 
used by the community yet, it is expected that this situation will change in the 
future, since it is the only interpretation that allows to explain systems composed 
by any number of star, and provide confidence intervals in the application of 
the models to real observations. The complete statistical formulation will be 
described in next section. 
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6.3. Complete statistical formulation 

In this section, we will explain the basics of the method of statistical modeling. 

For a complete descroption of the method, see lCerviho. Luridiana. Cerviho-Luridianal 
( 2005) . The general problem is the computation of the luminosity of an ensemble 
of stars. A luminosity distribution function (LDF: (p^i)) can be assumed, which 
describes the expected distribution of luminosity values in a generic ensemble. 
The integral of the LDF is normalized to 1: 

/•CO 

/ <p L (t)dl = 1. (13) 
Jo 

The integrated luminosity of an ensemble is traditionally obtained by means of 
the expression: 

/•CO 

L to t = N tot £<p L (£)d£. (14) 
Jo 

The bottom line of the statistical formulation is that the traditional ap- 
proach shown in Eq. ^] is conceptually wrong and operationally sterile. The 
crucial point here is the definition and interpretation of the LDF: (fh(£) is a prob- 
ability density function (PDF) from which the luminosities of individual sources 
are drawn; accordingly, the total luminosity of a cluster cannot be deterministi- 
cally computed, but its mean value M[ can be obtained as: 

M[ = (L tot ) = N tot / £<p L (£)d£. (15) 
Jo 

The integral on the right-hand side of this equation is the mean value of 
the LDF, fj,'i, so that: 

M[ = N totf i[. (16) 
In terms of the IMF and the isochrone, the LDF can be expressed as follows: 

^;t) = m (m)xf^)^ (17) 
V dm J 

where the time dependence of ipi, has been written explicitly. If we change the 
integration variable in Eq. EDfr° m & to m, the mean value of the LDF can be 
rewritten in terms of the isochrone and the IMF 2 : 



£(m;t) (p M (m) 



/ dl{m;t) \- 1 d£{m; t) 
\ dm J dm 



dm 



£(m; t) ipM(m) dm 



£i(t), 



(18) 



2 The isochrone is not monotonic, so that the integral limits of Eq. ll5l do not correspond to those 
of Eq. EU 
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where m and m up are the lower and upper mass limits respectively of the 
integration domain. 

Solving this integral is the main task of stellar population synthesis mod- 
eling. Note that, although Eq. ^] and the combination of Eq. [21 and Eq. |I] 
would lead to the same mathematical expression, the interpretation of what is 
computed differs drastically. In the statistical case, synthesis models compute 
the mean value of a probability distribution function (that can be split in com- 
ponents, or evolutionary phases), but, by virtue of the very probability concept, 
it is not needed that all the modeled phases should be present in an observed 
cluster. In the deterministic case, it is mandatory that all the considered phases 
are present (by construction), even if they correspond to an unphysical frac- 
tional number of stars (i.e fractional amounts). So, this statistical formulation 
includes by its very nature, the possible sampling effects in real stellar popula- 
tions. As an important point, it is necessary to note that sampling effects are 
not only a characteristic of the system under study, but also a characteristic of 
the observation (a narrow slit observation of a given system will include in the 
observation a lower number of stars than an observation with a broader slit). 

Additionally, Eq. ^] also reveals some of the problems and sources of un- 
certainty we have been talking about: 

1. The discontinuities in the derivative (d£(m; t) / dm)' 1 are directly related 
to fast evolutionary phases and discontinuities in the input tracks and 
isochrones. 

2. The fact that the models results come from an integral, /^ OT t ifM dm, 
instead of a sum of evolutionary phases, ^2iWi£i(t), implies the need of 
selecting actual representative evolutionary phases, and not the tabulated 
isochrone points as they are directly computed (especially if they are rep- 
resentative of a maximum luminosity, like the tip of the AGB). 

7. Conclusions 

In this review we have addressed the current status of evolutionary synthesis 
models. We have shown that there are several sources of uncertainties: 

• With respect to the input ingredients, the main source of uncertainty are 
interpolations. Such interpolations should ideally follow a physically-based 
scheme, but such scheme is not always possible. We have shown here 
the physical assumptions of current interpolation schemes for tracks and 
atmospheres. 

We have shown that the most reliable ages for model use and comparison 
are the MS turn-off ages of the used tracks. If differences appear at such 
ages, they reflect differences in the numerical methods used in each code, 
which must be further investigated. 

• With respect to the computations aspect of synthesis models, we have re- 
viewed the two main methods: isochrone synthesis and FCT. Since both 
methods produce different results, we concluded that the stellar popula- 
tion synthesis method are not a reliable tool yet. For the time being, 
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only those regions of the electromagnetic spectrum where both methods 
coincide should be relied on as realistic outputs. 

• Regarding the numerical computations performed by the models, we have 
shown that quantities like the SN rate (or the rate of any given evolution- 
ary phase) is a useful quantity to address both numerical and unphysical 
assumptions included in the models. 

• We have also reviewed additional uncertainties related to the incomplete- 
ness of the input ingredients, and we have outlined the issue of inclusion 
of stellar rotation in synthesis models. 

• Finally, we have discussed the use of synthesis models in realistic cases. 
We stress again that the results of the models depend also on the number 
of stars covered by the observation, and not only the number of stars in 
the system under study. We show how to address these problems with a 
statistical interpretation of the models' results. 

The main conclusion of this work is that the current literature about popu- 
lation synthesis does not provide enough information to address the uncertainties 
involved. In this sense, uncertainties can only be addressed if synthesis mod- 
els papers explain completely and carefully their input hypotheses, in particu- 
lar the interpolation schemes and its justification in physical or mathematical 
terms. An extensive documentation of codes (as Gary Ferland has done with his 
photoionization code) would highly improve the understanding and uncertainty 
evaluation in the synthesis method. 
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